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Abstract 

A general approach for deriving the expression of reference (density of) distribution functions, 
J-^, by statistical thermodynamics and the definition of local equilibrium conditions is illustrated. 
This procedure may be adopted for a system subject to an arbitrary number of thermodynamic 
forces. For concreteness, we analyze the case of a system submitted to three independent ther- 
modynamic forces and the local equilibrium corresponds to the configuration of minimum entropy 
production condition and the maximum entropy principle. In this limit case, we show that the 
derived expression of distribution function is more general than that one, which is currently used 
for fitting the numerical steady-state solution obtained by simulating the Ion Cyclotron Radiation 
Heating (ICRH) FAST-plasmas and for describing various scenarios of tokamak plasmas. Through 
kinetic theory, we fixed the free parameters linking them with the external power supplies. The 
singularity at low energy in the proposed distribution function is related to the intermittency in 
the turbulent plasma . As a matter of fact, this work is not restricted to, but applied to, tokamak 
plasmas. Tokamak-plasmas are taken as an example of close thermodynamic systems. An appli- 
cation to a simple model of fully ionized tokamak-plasmas submitted to an external Ohmic source 
is discussed. 

PACS Numbers: 52.25.Dg, 05.70.Ln, 05.20.Dd, 52.55.Fa. 

Keywords: Statistical Thermodynamics, Thermodynamics of Irreversible Processes, Kinetic The- 
ory, Tokamak-Plasmas. 
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I. INTRODUCTION 

Statistical thermodynamics constitutes a powerful tool for deriving the reference density 
distribution functions (DDF), J-"" (see, for instance, [T]). The DDF is obtained by starting 
from the Einstein formula linking the (density of) probability of a fluctuation, J-", with the 
entropy variation, AS, associated with the fluctuations from the state of equilibrium [2], [3] 
and |1] 

J^ = Afoexp[AS] (1) 

In Eq. ([I]), Boltzmann's constant is set equal to 1 and Ao is a normalization constant that 
ensures the sum of all probabilities equals one. Notice that, since AS is the change of entropy 
starting from equilibrium, we have AS* < 0. Prigogine generalized Eq. ([I|, which applies 
only to adiabatic or isothermal transformations, by introducing the entropy production due 
to fluctuations. Let us consider a system characterized by n variables Ai, ^2, ■ ■ ■ A„. The 
(local) equilibrium values are A^, Ag, ■ ■ ■ A°. Denoting by a^ = Ai — A'^ {i = 1 ■ ■ -n) the n 
deviations of the thermodynamic quantities from their equilibrium value, Prigogine proposed 
that the probability distribution of finding a state in which the values a, lie between a^ and 
tti + dai is given by 

J^ = Aroexp[-A^^] (2) 

where we have introduced the dimensionless [density of ) entropy production AjS [6] 

AjS =-cx>0 (3) 

n 

Here r, n and a denote the relaxation time, the number density and the entropy source 
strength (per unit volume), respectively. The negative sign in Eq. ^ is due to the fact 
that, during the processes, — A/S* < 0. Indeed, if —AjS were positive, the transformation 
ai — )■ a'i would be a spontaneous irreversible change and thus be incompatible with the 
assumption that the initial state is a (local) equilibrium state [S]. Using this notation, the 
thermodynamic forces X'^ and the thermodynamic fluxes J^ are 

X^=^ ; J. = ^ (4) 

aa^ dt 



and the expressions for djS and a read [7] 



diS = Y,X"dcik or ^ = ^J.X">0 (5) 






djS , 

adx. 



dt jQ 

where (ix is a (spatial) volume element, and the integration is over the entire space Q occupied 
by the plasma. Note that the probability density function ^ remains unaltered for flux- 
force transformations, X^ — )■ X'^ and J^ — !■ J^, leaving invariant the entropy production 

m 

K=l K=l 

The expressions for the reference distribution function, J-'°, can now be obtained by im- 
posing a set of minimum conditions able to identify, uniquely, the local equilibrium. These 
conditions should be established according to the particular physical situation that we are 
analyzing. It is important to stress that the set of restrictions, imposed by very general phys- 
ical principles, including scale invariance, gives automatically, for some range of parameters, 
a distribution function whose singularity can be interpreted in the terms of intermittency on 
turbulent plasma (see Appendix). This procedure is quite general and it can be applied to a 
system subject to an arbitrary number of thermodynamic forces. However for concreteness, 
we limit ourselves to study the particular case of a system submitted to only three thermo- 
dynamic forces. The local equilibrium condition is determined by imposing the minimum 
entropy production condition and the maximum entropy principle. 

The main objective of this paper is to derive the expression of J^^ for the case of a system 
subject to three independent thermodynamic forces by statistical thermodynamics. The link 
between the free parameters with the external power supply and the entropy source strength 
of the confined tokamak-plasma is performed by kinetic theory. This example of calculation 
illustrates how this procedure applies to the general case of systems subject to a larger 
number of thermodynamic forces. We mention that the present state of the art in tokamak 
modelling allows obtaining distribution functions from sophisticated numerical codes (e.g., 
NUBEAM, ASCOT, PENCIL for NBI, and TORIC, PION for ICRH). These numerical 
codes can actually be used for calculating numerically the particle distribution function for 
different power injection schemes and different levels of approximations in treating sources, 
collision operators and particle motion in the equilibrium fields. The scope of the present 



work is to provide a reasonably complete class of model distribution functions to be used 
in gyrokinetic or hybrid gyrokinetic-fluid simulation of nonlinear dynamics of drift Alfven 
waves in burning plasmas of fusion interest. The coupling of such gyrokinetic or hybrid 
gyrokinetic-fluid simulation with power deposition codes, like those mentioned above, is a 
challenging task, which is underway. In the meantime, sufficiently accurate description of 
model particle distribution functions are needed, which nowadays are typically chosen as 
"reasonable" model functions of the particle constants of motion (often just a Maxwellian, 
as in GTC and GYRO, or more generally an anisotropic Maxwellian or slowing down for 
the case of HMGC, M3D and NOVA). In this work, the constraints that non-equilibrium 
statistical mechanics impose on the adopted model distribution functions are derived and a 
class of model distribution function is proposed, whose usefulness is therefore, the readiness 
to be adopted in gyrokinetic or hybrid gyrokinetic-fiuid simulation of nonlinear dynamics of 
drift Alfven waves in burning plasmas of fusion interest. At the present time, none of the 
above-mentioned codes can be used in the same way, for the intrinsic difficulty of an actual 
integrated simulation. The advantage of the model distribution functions obtained in this 
work is therefore evident. As a matter of fact the work is not restricted to but applied to 
tokamak plasmas. The fundamental issue here is that the detailed application to tokamaks 
comes in only when the specific form of constants of motion in the equilibrium fields is 
adopted. Before that point, there is a construction of the particle distribution function 
out of the equilibrium, as we expect that to be in the case of a tokamak. The final part 
of the work is addressed to study the inverse problem: given J-"", we are looking for the 
Fokker-Planck equation admitting J-"" as the steady-state solution. 

The paper is organized as follows. In Section (III]) we derive the general expression of the 
reference density of distribution probability, J-"", by a purely thermodynamic approach. The 



parameters, entering in the expression of J-''', are determined by kinetic theory in Section III 



by adopting a model for the tokamak-plasma and the external sources. Section [TV] addresses 
the following questions : 



For collisional tokamak-plasmas, how much is the deviation of the reference DDF, J^^ , 
from the Maxwellian ? 

Does this deviation coincide with the one estimated by the neoclassical theory (see, for 
example, Ref. [12]) ? 



These are crucial questions because, when a new approach is proposed, one of the main 
tasks is to show that the new approach is also able to recover the results, which are well 
estabhshed in literature. In our case, we should ensure that the expression that we found for 
the reference DDF, J-"", coincides exactly with the one predicted by the neoclassical theory 
for collisional tokamak-plasmas. As we shall see in this question, the answer to this question 
is positive. Some concluding remarks can be found in the Section [V] 

II. THERMODYNAMIC DERIVATION OF THE DISTRIBUTION FUNCTION. 
GENERAL CONSIDERATIONS. 

In this section, we derive the form of J-"", the reference density of probability function, by 
a pure thermodynamic approach. We limit ourselves to study a system subject to three 
thermodynamic forces. As usual, the gyrokinetic (GK) theory makes often use of an initial 
distribution function of guiding centers. In the GK simulations, as well as in the GK theory, 
this initial distribution function is usually taken as a reference DDF if it depends only on 
the invariants of motion and it evolves slowly from the local equilibrium state i.e., in such 
a way that the guiding centers remain confined for sufficiently long time. After a short 
transition time, the state of the plasma remains close to the reference state, which results 
to be a small deviation of the (local) equilibrium state. The expression of the coefficients 
of the J-'° will be determined in the next section by kinetic theory. The reference DDF is 
obtained by perturbing the local equilibrium state. To determine the reference DDF we 
shall follow the following strategy. First we consider open thermodynamic systems obeying 
to Prigogine statistical thermodynamics. Successively, we define the local equilibrium state 
by adopting a minimal number of hypotheses. Finally, we link the Prigogine probability 
distribution function with particles DDF. The expression of the local equilibrium state is 
fixed by adopting a set of minimal number of ad hoc assumptions. We restrict ourselves by 
analyzing the case of an axisymmetric magnetically confined plasma. After having performed 
the guiding center transformation, the necessary variables for describing the system reduce to 
four independent variables [6]. In addition, we consider the case where the system is subject 
to three thermodynamic forces, associated to the three fiuctuations ai, 02 and a^ with a^ 
identified as w, the particle kinetic energy per unit mass. Hence, w = {v'n + v'j_)/2 where 
f II denotes the parallel component of particle's velocity, which may actually be parallel or 



antiparallel to the magnetic field and v± is the absolute value of the perpendicular velocity. 
Now we should define the local equilibrium state. We assume that (see also the footnote 

m) 

The equilibrium state corresponds to the values ai = and a2 = for which the entropy 
production tends to reach an extreme. 

Hence, ai and 0:2 are two Prigogine's fluctuations [for their exact definition see Ref. j5] and 
[13] . For fully ionized coUisional plasma, these fluctuations are linked to the thermody- 



namic forces by Eq. ([59j) (see the Section IV)]. At this state, the deviation of AjS from its 
equilibrium value is given by [3] and [13] 

— AjS = go{w) — - y^ gij{w)aiaj + h.o.t. with (7) 

go{w) = -AjS \ai=a2=0 

where h.o.t. stands for higher order terms. Therefore, the general DDF becomes J-"*^ when 
the expression of entropy production is given by Eq. ([T]). Coefficients gij denote the second 
derivatives of AjS with respect to the a, variables. With these coefficients we may form a 
positive definite matrix. This matrix can be diagonalized by obtaining 

2 
- AjS = go{w) - J2 ^^H6,J^■ + h.o.t. (8) 

where Ci{w)Sij is a positive definite matrix and 5ij denotes Kronecker's delta. Eq. ([s]) al- 
lows describing the entire process in terms of two independent processes linked to the two 
independent fiuctuations ^1 and ^2- As mentioned, J-'° is determined by four variables. In 
order to simplify, as much as possible, the expression of the reference density of probability 
function, we shall now introduce a set of convenient variables. We introduce two variables, 
Ci and (2, of the type, degrees of advancement (for a rigorous definition, see any standard 
textbook such as, for example, Ref. [13]. See also the footnote [32]), and the remaining 
variable, (3, is chosen such that the fiuctuations .^1 and ^2 depend on the variables (i and 
(2, but not on variable (3. The reference density of probability function will be then deter- 
mined by the following set of coordinates: {w, (1X2,(3)- We recall that the remaining two 
coordinates: the toroidal angle, 0, and the gyro-phase angle, $, are unnecessary due to the 
axisymmetric configuration and the guiding center transformation, (i and (2 are associated 
to the two independent fiuctuations ^1 and 1^2 (therefore, they are also associated to the two 



independent thermodynamic forces). We may expand the fluctuations ^i and ^2 in Taylor's 
series of variables (Ci, C2) by obtaining 

In Eqs. ([9]) we have taken into account that Ci and C2 are two degrees of advancement 
variables and that ^1 and ^2 are two independent fluctuations (see the footnote [33]). Hence, 
the reference distribution function ^ takes the form 

dJ-o = V{w) exp[-ci(w)(Ci - C?)'] exp[-C2{w){C2 - C2f]dT with (10) 

V{w) = exp[5fo(w)] = exp[-A/S']g,=^2=o and q(w) = Ci{w)[-^] ^ ^ (^ = 1,2) 

The phase space volume element dT = d'^dw is linked to the new volume element dV = 
dwdC,id(^2dCzd(j)d^ by 

dT=\J\df (11) 

where | JT | is the Jacobian between dV and dV . To proceed further in our analysis, we should 
complete our definition of equilibrium state. For this, we may think to u; as a random variable 
satisfying two restrictions 

E[w] ~ T/m > (12) 

E[ln(u;)] ~ ^{T/m) 

where T is temperature and m denotes particle's mass. E[ ] and "^{x) are the expectation 
operation and the digamma function, respectively. The equilibrium state is fully determined 
by assuming the validity of the following Principle of Maximum Entropy. 
In the phase space T, at the externalizing values ai = 0:2 = (or, which is the same, at 
^1 = ^2 = 0), the system tends to evolve towards maximal entropy configurations under the 



constraints given by Eqs (12) 



Under this assumption, the random variable w is modeled by a gamma distribution 



V{w/e) = KQY' exp[-w/&] (13) 



where we have introduced the scale parameter B and the shape parameter 7. Conditions (12) 



fix the parameters of the Gamma distribution (see, for example, Ref. [15]). The motivation 



for the choice of the restrictions given by Eqs(12), as well as the special mathematical 



properties of the resulting PDF, are given in the Appendix. Hence, the maximum entropy 



principle allows writing the reference state, Eq. (10), as 



-s / in\ 7^1 

d^ = ^o[l) exp[-^/e] exp[-ci(^/e)(Ci - C?)'] exp[-C2(u;/e)(C2 - (2°)'] I J I ^r (14) 

If we interpret our reference DDF as a time and ensemble average of the physical DDF 
describing turbulent plasma, then the singularity at w = for < 7 < 1 can be related to 
the intermittency (see Appendix). Notice that at the point with coordinates (0:1,02, as) = 
[0,0,(7 ~ 1)0]) the system satisfies the principle of maximum entropy and the entropy 
production reaches its extreme value. 

Let us now consider the coefficients Ci{w/Q) and 02(^/6). As an example of a closed 
thermodynamic system (i.e., the system can exchange energy but not matter with an outside 
system), we analyze collisional tokamak-plasmas where the magnetic configuration can be 
brought into the form (see, for instance, Ref. [12] ) 

B = -^ViJxV(j) + F{ij)V4> (15) 

zvr 

The particle energy is f = mw + qAe ~ mw, with m, q and Ae denoting the mass particle, 
the charge particle and the electric potential, respectively. Moreover, ip indicates the poloidal 
magnetic flux and F is the characteristic of axisymmetric toroidal field depending on the 
surface function tp. We identify (i, (2 and Cs with the toroidal angular moment, P^, the 
variable. A, and ip, respectively. These quantities are defined as (for a rigorous definition, 
see any standard textbook such as, for example, |16j ) 

Bo Fv\i , , 

"Oc I -D I 

sin^ dp II v\ 

A = = — With /i 



2\ B\ w "^ 2 I 5 I 

Here Qqc is the cyclotron frequency associated with the magnetic field along the magnetic 
axis. Bo- \ B \ and 6p denote the magnetic field intensity and the pitch angle, respectively. 
We expand these coefficients in terms of w. Considering that w/Q behaves as an entropy, 
and that (P^ — P<^o)^ behaves as f j?, for coherence with Eqs (1 ) and (2), the expansions should 
be truncated at the zero order for ci and at the first order for C2 (Ci and C2 are evaluated 



at the zero order). Hence 



ci(w/0) = c!°^ + /i.o.t. 



{l)W , 

Cn ' + c\' — + h.o.t. 
^ ^ 



(17) 



The expressions for Ci and C2 can be brought in the more convenient form 

1 n2 






SO) 



AAo 



+ 



w 



const. 
> 



AAoVAAi 0. 
Finally, the expression for the (density of) probability distribution function J-"*^ reads 



(18) 



J''=^fo{'^y~eM-w/&]exp 



P(t> ~ P<t> 



where A/q ensures normalization to unity 



exp 



j^"dr = 1 



AAq w\ (A-Ao)^ 
AAi q) AAo 



J\ (19) 



(20) 



.(0) 



'Q, Jn 

Observe that the presence of the free parameter €"2' is crucial. Indeed, as we shall show in 
in particular, see Eq. (81)], the absence of ^ precludes the possibility of 



the Section 



IV 



identifying the DDF, given by Eq. (19), with the one estimated by the neoclassical theory 



for collisional tokamak-plasmas (see, for example, Ref. p2j)- In addition, it allows describing 
more complex physical scenarios such as, for example, the modified hi Maxwelian distribution 
function [17J. Last and not least, in some physical circumstances, the presence of C2 is essen- 
tial to ensure the normalization of the DDF. Thermodynamics has been able to determine the 
shape of the DDF, but it is unable to fix the seven parameters 0, 7, P^o? Aq, AP^, AAq, AAi. 
These coefficients have to be calculated in the usual way by kinetic theory. Figs ([T])-(|2]), ([s])- 
(|4]), and ([5])-(|6| illustrate three plots of Eq. (19) (estimated for unit values of the Jacobian 
and the normalization coefficient) corresponding to the values w = EQ, 'j = 1 + E (with E 
denoting the Euler number), P^ = P^q and A = Aq, respectively. 



III. DETERMINATION OF THE PARAMETERS 

This section is devoted to the determination of the value of the parameters appearing in 
Eq. (19) for J-"^. As an example of application, we analyze tokamak-plasmas in collisional 



transport regimes. To accomplish this task we should 
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FIG. 1: Distribution function, Eq. (19), (per 



FIG. 2: Contour plot of Fig. [T} 



unit values of the Jacobian and the normalization 
factor) computed aX ^ = 1 + E [E indicates the 
Euler number), w = EQ, AP^ = 22.360, AAq = 
50.00 and AAi = 30.2031. 

• Adopt a model for tokamak-plasmas; 

• Model the source terms. 

Calculations should be performed by a kinetic approach. For definiteness, the particle 
density Ua and the average velocity Uq, of each species a (with a = e for electrons and 
a = i ioT ions) satisfy the equations 



n„ X 



rfvj^""(v,x) 



(21) 



n„ X u X 



/ t/v vj^"°(v,x) with a = {e,i) (22) 

Jv 

where V is the velocity- volume in the phase-space. In addition, we should construct repre- 



sentation ( 19 ) in such a way that temperature T^ does coincide with the exact value of the 



temperature of each species. This implies 

1 



na(x)Ta(x) = -rUc 



dv\\r-Ua r J'°"(v,x) 



with 



a 



[e,i) 



(23) 
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FIG. 3: Distribution function, Eq. ([19]), (per FIG. 4: Contour plot of Fig. [3} 



unit values of the Jacobian and the normalization 
factor) computed at 7 = 1+E, AP,^ = 22.360 and 
A = Ao. 
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FIG. 5: Distribution function, Eq. (19), (with FIG. 6: Contour plot of Fig. [5} 



the Jacobian and the normalization factor set to 
1) computed at ^ = 1 + E, P^ = P^q, AAq = 
11.111 and AAi = 50.00. 
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where m^ is the mass particle of specie a. To the Eqs (21 )-(23 ), we should add the equations 
warranting the validity of the energy and entropy balances. To this end, we consider three 
kind of external energetic sources: the total Ohmic heat, ^S-q^-, supplied from outside, the 
energy A£^^^j, supplied by neutral beam injection (NBI), and the energy of the minority 
population heated by ion cyclotron resonant heating (ICRH) ^£f(j^jj. Here, we shall not 
deal with burning fusion plasmas and the loss due to the Bremsstrahlung effect is neglected. 
At the end, the energy balance equation reads 

im„ / rfr I V - u, p ^"0(v, x) = A^S,. + A£°^, + ^E'^crh ' ^^ (24) 

where A£^£ is the energy loss and we have neglected the electric contribution with respect to 
the kinetic energy. Let us now consider the entropy balance. The total entropy, s, is defined 
as [6] 



nr,[^]Sr 



[ rfv^°(v,x)ln('^^0(v)) = - / rfv^O(v,x)ln^O(v) (25) 



Here h is the Planck constant, e is the (positive) magnitude of the electric charge and m,a 
the mass of the species a. This equation is, however, inconvenient for setting the values of 
the parameters, because it involves the total entropy density that, at this stage, is not a 
known quantity. It is more useful to study separately the two entropy contributions, (IeS'^ 
and diS"", of the total entropy S"": 

(i5° = rfs5" + diS" with 5" = /" rfx n«(x)s«(x) (26) 

Jn 

(IeS"' represents the amounts of entropy crossing the boundaries. The total entropy flux 

of species a is given by two contributions: the convective entropy flux, n^SaU" and the 

conductive entropy flux, Jg. On the contrary, djS"" represents the source entropy, due to 

internal processes, which remains within the system. In terms of the conductive entropy 



flux, J^, and the entropy source strength, a", Eq. (26) takes the form 
dS" 



- f dA- [n„s„u" + Jg(x)] + f dx (t"(x) (27) 

is Jq 



dt 
where S denotes the boundary. According to the kinetic theory we have [7j 



Jg(x) = - / dv [v-u"(x)]J^"°(v,x)ln( -J-"°(v,x) 

= - f d^j[^j- u"(x)]J^"°(v,x)lnJ^°°(v,x) (28) 

Jv 
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and 



a" 






(29) 



where /C"'^ denotes the colhsional operator of species a due to (3 and we have taken into 
account Eq. (|3|. From Eq. (28), we obtain 



- / rfv [v - u"(x)]^"°(v, x) ln^"°(v, x) = hje, - Jso.. ' ^e^,, - Js.cuh) (30) 

Here J^:^ indicates the energy loss flux. J^q,^ , J^-jvsj ^^^ ^£icrh denote the Ohmic energy 
flux, the NBI energy flux and the ICRH energy flux, respectively. T^ is the temperature of 



species a. Eqs (24), (29) and (30), together with the definitions given by Eqs (21 ) and (22) 



allow determining the seven parameters 6, 7, P^o? -^O) AP^, AAq, AAi. More particularly, the 
values of these parameters may be fixed by adopting the following strategy. The coefficients 



7, AP(i, AAo and AAi are provided by Eq. (29), whereas Eq. (24) sets up the calculations for 



the determination of parameter B. P^o and Aq are determined by expansions ([9]). The plasma 



temperature profile can be obtained by Eq. (30). This solution should be in agreement with 



the definition given by Eq. (23). Generally, the latter calculation is very complex, but it can 



be strongly simplified by adopting the following procedure. In a first phase, the plasma is 
heated to thermonuclear conditions by the externally supplied power. In our simplified case, 
the applied sources are the sum of Af^o/i., ^£nbi and ASjcrh- When the thermonuclear 
conditions have been achieved, the plasma temperature can be sustained if, in the overall 
power balance, the power loss, Vl, is balanced by the external supplied Ohmic (Voh.), NBI 
(Vnbi), and ICRH (Vicrh) powers as well as the alpha-power heating {V alpha) i-e.. 



Valpha + 'Poh. + 'PnBI + 'PiCRH — 'Pl 

Jg(x) = 



(31) 



Of course, Eq. (24) and Eq. (30) remain valid also for burning plasma, which can be used 
after having specified a proper model for estimating the source Vaipha- However, as al- 
ready mentioned, the loss due to the process of producing the radiation (Bremsstrahlung) 
is neglected and Vaipha is not treated in this work. Since the quantities entering in this 



equation depend on the plasma temperature, Eq. (31) may be solved for estimating the 



average temperature. However, as previously mentioned, calculations can be performed only 
after having modeled the external energy sources and estimated the entropy source strength 
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a" by modeling the magnetically confined plasmas. Below, we provide some examples of 
modeling. 

The Ohmic Heating Model 

The expression of the Ohmic heating density can be found in many reference books. For 



the magnetic configuration (49), we have (see, for example, Refs [16] and [T8] ) 



p ^ -2 = 9 [ ^ ^° 1^ \ ('X'y\ 

I Oh. m ''[47ri?o(l + r/i?ocos^)j g(0)[g(a) - l/2g(0)] ^ '' 

where rj is the resistivity of the plasma and j is the current density, a and rjs denote the 
minor radius of the tokamak and the Spitzer resistivity, respectively. 

The NBI Source Model 

Once the neutral beam enters in the plasma, the neutral particles will be ionized. Their 
energy passes to the particles of the plasmas causing heating of both electrons and ions. The 
energy supplied by neutral beam injection may be modeled thinking in terms of a pencil 
beam source, S^bi, acting for a period of time r 



Snbi = ribrSiy - v^) (34) 

where nAr^/ is the birth rate per unit volume and v^ is the beam ion velocity. Hence, the 
NBI energy is 

^^NBi - -TzmbfibT I Vfe 1^ fip (35) 

where m^ and Vtp are the the mass of the beam ion and the plasma volume, respectively. 

The ICRH Source Model 

Let us consider the case of a low concentration of ions ^He colliding with a thermal 
background plasma, composed by Deuterium and electrons. The ^He minority is about 
2% — 3% of the density of the background plasma and it is heated by ion cyclotron resonant 
heating (ICRH). In the velocity space, the long term evolution of the distribution function 
for the high frequency heated ions, is governed by [20], (see the footnote [51] ) 

^^'^ = -V ■ S'"(J^'") + P(J^'") (36) 
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where 

S™(J^"^) = S;^(J^™) + J2 Sri-^"") (37) 

a=e,i 

Suffices m and a [with m = ^He and a = (e, i)] distinguish the minority population and 



the species of the background plasma, respectively. The ffist term in Eq. (37) describes the 
quasi-linear diffusion due to the resonant wave particle interactions and the second term 
is due to the coUisional operator. P(J-'™) takes into account other auxiliary sources; in 
our analysis, we shall put P(J^"^) = 0. The gradient operator, V, is defined as the row 
vector V = {dv^,dvy,dv^) in the velocity space, whilst V ■ A is the matrix multiplication 
between the gradient vector and the matrix A. In this approach, the ICRH is modeled by 
the divergence of the quasi-linear flux S|^(J-'™), which is proportional to the energy density 
of the propagating wave 

ASicRH ^ -\m„.r [ dT \ ^r - u^ \' V ■ S™ (^"^) (38) 

^ J{nv} 

where r, m^ and u^ are the small laps of time during which the minority population is 
heated by ICRH, the mass and the average velocity of the minority, respectively. S™" and 
the simplified quasi- linear term, S™, can be written as 



gma^jr^) = _[V • (d""^'V"^)]^ + D™"(i)j^™ (39) 

with T denoting the transpose operation. We introduce the following dimensionless coordi- 
nates, w and A 



w = -g— where Vtha = \ — - « = (e, i) (40) 

X = BnX 
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with w and A defined in Eqs (16)]. In these coordinates, the matrices D^ ' ', D^y' ', D 



mo(2) 



™(2) -f^ma{l) 



and D^ can be cast into the form 



= me(2) 



■Qme(l) 



■j-jmi(l) 




= mi(2) 



T^me{2) _|_ 2V^ dP, 



+ 



me(2) \ 



Vthe 

me{2) 



dw 



1-2\B\X p, 
2|B|A«) cAA 

1-2|B|A jjmi(2) 
2\B\\w cAA 



D 



mi{2) 
cww 





D 



mi{2) 
cAA 



(41) 



with I i? 1 = 1 i? I /i?o, and 



D 



m(2) 



T-)m(2) T-)m(2) 
-^VFAw -^VKAA 



/ 



D 



m(l) 



_2^p)rn{2) 2v^ '^^I^l^l i ^_ 

2 jjM2) , 2v^ '^-Pw-aJ I 1 

\ VtheVli W^"^ Vthe dw Vthe ' 



(42) 



l-2|i?|A jyn{2) ^_ \ 2\{\-2\B\\) dP'^^jX 
2\B\\w ^^^ Vthe Y \B\w dX 



l l'2\B\X p,m(2) ^ / 2A(1-2|J3|A) rfD;;;'^^ 



dX 



J 



Notice that in Eqs (41) and (42), extra terms appear in the components of the drift vectors. 



However, the nature of these additional terms is purely geometrical and they come from the 



conversion of the balance equation, Eq. (36), into the standard Fokker-Planck equation (see, 
for instance, Refs [23] or [23]]. Hence, the expressions of the coUisional electron and ion 
coefficients can be brought into the form 

'erf (v^) — V^ erf( 



jjme{2) 



D 



me(2) 
cAA 



2VtheVw 
pme 



W 



(43) 



((2-i)erf(v^) + 



erfVwl 



7-)me(l) 
cw 

-pma 

-^Dma ~ 



T"^^ rrii [erf (a/^) — \/ili erf'(\/M))] 



^the me 

4V2^nae^ZlZl\nAr 
3 m'L 



where 



w 



Ar 



3(Tm + TajXpr, 

AZjyjiZjQ^e 



^m^ oA^m ~T ^a) 



^'K Zj,p[iZjfy_e \TiyYiJ-m ~T rioil Qi^ 



1/2 
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and 



jjmi{2) 



D 



mi{2) 



2VtheVw 



erf (vKti)) — y/Kw erf ( 



KW) 

1 \ „, ^^. erf'f 
— i< 



;/ 



erf(v Kw) + 



KW 



jjmiil) 



r™ (erf(-\A«£;) — \/~k^ erf'(-\//«i))) 



where 



V 



the 



W 



K = 



-'the 



Hhi 



(44) 



respectively. Notice that, the quantities in Eqs (43), are expressed in the CGS Gaussian 



units, erf(x) and erf (x) denote the error function and the derivative of the error function 
with respect to its argument, respectively. The expressions of the diffusion coefficients for 
the resonant wave particle interactions are 



D 

D 
D 



m(2) 
Www 

m(2) 
WwX " 

m(2) 
WW ~ 



2\ B\ X D 



±± 



(45) 



D 



m(2) 
WXw 



(2 I 5 I A)(l-2 \B\\) D] 



;i - 2 I fi I A) D^^ 



DT, = D^Jl{Usj2 \B\Xw) 



±± — U "p ' 
Vthel^l. 



Vtr 



Qr 



Zui^Bq 



Tn,fYiC 



with 

P^^^ ~ 10 ^ 80 W/cm^ 



^0 



plin 
abs 



AmmUm /o°° x^Jp{imx) exp (-x^) dx 



Here, Jp{x) indicates the Bessel functions of the first kind, e and c are the absolute value of 



the charge of the electron and the speed of light, respectively. From Eqs (41 )-(45 ) we obtain 



the expressions of the total matrix diffusion coefficients and the total drift vector coefficients 

I T-)m(2) r^m{ 



D 



m(2) 



= me(2) =mi(2) =m(2) 



D' 



D 



^wX ^XX 



(2)\ 
wX 

m(2) 



(46) 



J3m(l) ^ J3me(l) _^ J3mi(l) _^ j-)m{l) 



'W 



D 
D 



m(l) 
w 

m(l) 



In the zero orbit width limit, the steady state DDF of the Fokker-Planck equation, with the 



diffusion and drift coefficients given by Eqs (46), may be fitted by the profile 



J^"^OC 



e 



exp 






w 



exp[--j6{X-Xo) \J\ 



(47) 
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fC''<x yo^) 




(i.2Si 



^ m(2) 

FIG. 7: Dimensionless diffusion component D^w versus the normalized, dimensionless variables 



w and A. The dimensionless diffusion coefficients are defined as D 



m(2) _ T-,m(2) ,-pme 



«J 



n ^ 'at /'pme 



V 



E/^l^^x i<y^ 1 I 




FIG. 8: Dimensionless diffusion component D ^ versus w and A. The matrix of the diffusion 



coefficients is symmetric: D, 



m{2) 
wX 



D 



m(2) 
Xw 



with 5{\ — Ao) indicating the Dirac function [25]. Figs ([T]), (tsl) and dol), report on the 



graphics of the dimensionless elements of the total diffusion matrix, D 



m(2) Am(2) 
«'«' ' ^wX 



D 



m(2) 
Xw 



'XX 



and D™ [see the first equation in Eqs (46)], versus w and A, respectively. The dimensionless 



.m(2) 



diffusion coefficients are defined as Z)™^ = D™" Vthe/'^"^'^- Figs (lO)-(ll) report on the 



m(2) 



dimensionless components of the total drift coefficients [see the second equation in Eqs (46) 
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Ef^^^^iX ]il^') 




FIG. 9: Dimensionless diffusion component -D^ versus w and A 



D^'"« l(J*) 




9:30 



FIG. 10: Plot of the first dimensionless component of the total drift vector, D^ , versus w and 

defined as D^ = D^ v^f^^/V^^, against w and A. These pictures have been obtained 
by considering only the fundamental cyclotron heating, p = 0, (minority heating) and by 
setting K,± ~ p^^ (with p^^ denoting the Larmor gyroradius). Moreover, as an example 
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^ mil) 

FIG. 11: Plot of the second dimensionless component of the total drift vector, D^ , versus w 
and "A. 

of calculation, we chose the following values of the tokamak parameters: Bq = 3.45 T, 
Rq = 2,96 m and a = 1, 25 m. Moreover, P^^" = 50 W/cm^ and n^ = 2.5% rij, respectively 
and, as an example of calculation, the fields have been estimated at the center of the tokamak. 
Notice that in our case, the charge number and the ion mass are Z^ = 1 and rrii = m£, = 2mH 
for the background, and Zm = Z SHe = 2 and rrim = m 3He — ^Tnu (with inn indicating the 
mass of proton) for the minority, respectively. 

ThePower Loss 

Neglecting the Bremsstrahlung loss, the rate of energy loss is mainly due to the thermal 
conduction losses. Its expression can be found in many reference books [see, for example, 
Refs [16] and [18]]). We have 



Vl = snp- 



TlT 



(48) 



where n and T are the average density and the average temperature, respectively, te is the 
energy confinement time. 

Modeling Tokamak-Plasmas 
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The entropy source strength a" can be estimated by modehng the tokamak-plasma. For 
instance, let us consider fully ionized tokamak-plasmas, defined as a collection of magnet- 
ically confined electrons and positively charged ions. In the local triad (6^,69, e,^), the 
magnetic field, in the standard high aspect ratio, low beta (the plasma pressure normalized 
to the magnetic field strength), circular tokamak equilibrium model, reads (see, for example, 
Ref. [12j) (see the footnote [35]) 

^ r Bo . 

g(r) i?o 1 + {r/Ro)cos9^ ^ ^ 

Here Bq is a constant having the dimension of a magnetic field intensity, and q{r) and Rq 

are the safety factor (in Ref. [I2] the reader can find an exact definition of this tokamak 

parameter) and the major radius of the tokamak, respectively. In the magnetic configuration. 



given by Eq. ( 49 ) , we have 



F = BoRo ; V(r) = 2nBo / — - dr (50) 

Jo qij) 

The dimensionless entropy production of species a, A/S*", is derived under the sole assump- 
tion that the state of the quiescent plasma is not too far from the reference local Maxwellian. 
In the local dynamical triad |26], A/S*" can be brought into the form (see Refs [12] and |27j ) 

A,^^ = g!ii(,« - 4'') + q^^gf^ + -gf^) + gjl^.W - gf^) + q^\gf^ + gf^) 

+ <lt'<'' + €^'3r + %T9f^ (51) 

^i^ - %js Uii + fi'ii ) + Q\\b {9\\ +g\\ ) + Q\\b 9\\ + Qpd 9p 

Here g" [with r =(p, ||, A)] denote the Hermitian moments of the distribution functions 
and g? , g? are the dimensionless source terms. Index n takes the values n = (1,3,5). 
In the linear Onsager region, and up to the second order of the drift parameter e, it can be 



shown that Eqs (51) simplify to [12] 

A,^^ = ai,(^« - -gf^r + mgf^ + -gf'^r + m-gf'^r + ^a^igp - grw + % 

+ 2711 (<^ - -gf^Ygf^ + il\{gt^ + gt^)gt^ 

+ -A9'^'y + -li9f'r-2a^9r9f' (52) 

AjS^ = f^{gf^ + sl'^r + m^^r + 25^(^f + -gl'Y9T^ + R\{gf^f 

where coefficients 5",., <5r-, ^^r indicate the dimensionless component of the electronic con- 
ductivity, the thermoelectric coefficient and the electric {a = e) or ion {a = i) thermal 
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conductivity, respectively. Moreover, 7|j, (5? and e? are the parallel transport coefficients in 



21M approximation. As shown in Eq. (29), the entropy production is closely associated with 



the collision term. Notice that Eqs. (52) have been derived by using the following linearized 
Landau collisional term 
2TcelellnA 



/Cf 



■-a'='l3 ■ 



rrir 






+ /f (V2) 



d Id 

-/o"(vi)-/o"(vi) 



c/3 



/o'(v2)-/r(vi: 



d 



/o'(v2)) (53) 



Moreover, Grs{g) and In A denote the Landau tensor and the Coulomb logarithm, respec- 
tively 



Grsis) 



g^Srs - QrQs 



In A 



3/2(Te + 7^.)Ag 
Ze2 



A 



f 4:71 Ze^ in fTe + UiTi 



D 



l 



TeT,(l + Z) 



-1/2 



(54) 



Here Z is the charge number of ions. In addition, in Eq. (53), the distribution function is 
expanded in powers of the drift parameter e 



r(v) = /o"(v) + e/r(v) + 



with 



/o"(v) 



"l27rT„ 



771, 



exp 



V 2r„ 



a 2 



(55) 



By summarizing, the expression of the entropy source strength, Eq. (52), may be put on the 



left-hand-side of Eq. ( 29 ) provide that the integral appearing on the right-hand-side of this 



equation is evaluated by means of the linearized Landau operator Eq. (53). 



IV. ESTIMATION OF THE SEVEN PARAMETERS APPEARING IN THE EX- 
PRESSION OF J^ FOR THE SIMPLE MODEL OF FULLY IONIZED TOKAMAK- 

PLASMAS 



As an example of calculation we estimate now the seven parameters of J-"" when the entropy 



production is given by Eq. (52). In particular, our purpose is to link the expressions of the 



coefficients Qj,^, appearing in Eq. ([7|, in terms of the transport coefficients. In addition we 
shall prove that the expression that we found for the reference DDF, J-'°, coincides exactly 
with the one predicted by the neoclassical theory for collisional tokamak-plasmas. 
From Eqs ([2]), ([T]) and (13), we have 

A n I i\i ('^\ W \ 2 1 2 /i-^\ 

A/S = -(7 - 1) ln(^Q j + e + 2^^^"^ ^ 2^^^"2 + 9\20L\a2 (56) 
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By expanding the previous expression around the reference value w = wq we obtain, up to 
the second order 



AjS 



I ^^, (Wq\ Wo 



+ 



-(7 



^ 1 ll 



{w - Wo) + (7 - l)7r^iuj - woY 
2w^ 



1 2 1 2 

-giia^ + -fi'22a2 + 5'i2aia2 + h.o.t. 



(57) 



Up to a normahzation constant, we have that the distribution function Eq. (19) is approx- 



imated by a Gaussian density distribution function (in the variable w) by setting to zero 
the coefficient of the linear term [i.e., —(7 — l)/wo + 1/6 = 0]. The global optimality 
conditions is obtained by imposing that also the sum of the constant terms vanishes [i.e., 
— (7 — 1) ln{wo/Q) + Wo/Q = 0]. These two requirements are simultaneously satisfied only if 

7 = 1 + ^ (58) 

Wo = (7 - 1)6 = EQ 



Hence, solutions (58) ensure not only a local approximation, valid up to the second order, but 



also a good global approximation. In this sense the values of Wq and 7, provided by Eq. (58), 



are optimal. From Eq. (H) and Eq. (57) we obtain the expressions of the thermodynamic 
forces 



^1 _ dAjS 

A = -^ = gnai + gi2a2 

^2 _ dAjS 

^ = -T. = gi2ai + fi'22a2 



da-) 



(59) 



X 



3 _ dAiS 



[w - Wo) 



dw ^92 
By solving the previous system of equations with respect to ai, 0:2 and {w — Wq), we find 

Oil = g22X^ - gi2X'^ 

012 = -gi2X^ + giiX'^ 

{w - Wo) = EQ'^X^ where 



(60) 



9j^ 



9jK 



[with (j,«:) = (1,2)] 



g = giig22 - g^ 



1 
g 



g gug22 - gh 

Hence, in terms of the thermodynamic forces, the electron and ion entropy source strength 
read, respectively 



AjS'^ = \EQlxf + ]^gl,Xl' + ]^gl,xf - gl.XlXJ + h.o.t. 
AjS^ = \Eejxf + \g\,xf + h.o.t. 



(61) 
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By diagonalizing, the first expression of Eq. (61) can be brouglit into tlie form 

1 



AjS" 



/"1,2 



-mlxf + ^^^ee' + l^2if 



with 



{hi + 922) ± V (^11 - 922^ + 4^ 



12 



ie = -j;9i2X^ - —{gu - 922 - h)X^ 



q2 



92 



L = -J912X^ + —(^11 - ^22 + h)X^ 

h = \ {gu - 922? + 4^?2 



?1,2 



\/4^?2 + (^11 - fi'22)(^ll - ^22 ± h) 



(62) 



(63) 



I ^12 I v^ 

Let us now evaluate the coefficients ^^^ and Gq, for the case of fully ionized tokamak-plasmas 
described, for example, in Ref. p2]. In the banana regime, without considering the classical 
contributions, the expressions of the electron and ion entropy source strength take the form 



AjS" 



.22 



.32 



6f zr + qzr + a\\z'r + ii\z\zi + i^^zizi + 2511^,^^,^ 



~AZ}'' ^~^e.zf ^2b\Z]zJ 



"iri 



(64) 



where 



^1! 



e(5) 



e(3) , -e(3) 



^e — 9\\ + 9\ 



zl^g 



(1) ^e(l) 

^11 



^ = 9 



i(3) , -j(3) 



(65) 



As mentioned in the introduction [see Eq. p])], there exists the following (linear) thermody- 
namic covariant transformation (TCT) (see Eq. ([6]) and Ref. [8]) 



XI 



-e(5) 
^11 



\^2 _ 72 _ „e(3) -e(3) 



such that 



X} 



i(3) , -i(3) 



Xj — aiZ^ 



Zl = gf + gli 



(66) 



AiS^ = AjS"' 
AjS' = AiS'' 



(67) 



The transport coefficients g°j and the transformation-coefficients aa, h^ [with a = (e, i)] and 



25 



Ce, are easily obtained from Eqs (61), (64) and Eqs (67). We find 
. _ 271, . , 



2aii 



2./a\\ 



\fEd\\Q, 



EQ, 



2__ 2 2 

9u = — ('^i^ll - "J) ; 922 = — (ef^ll - 7||) ; 9i2 = — ("||7|| - 5p||) (68) 



ail 



o"ii 



cm 



2Si 



2jKi 



E^Qi 



EQi 



) yii 



9u 



(eU^ -6?) 



By inserting the values of the collision matrix elements (see, for example, Ref. [6]), we find 



0.1168 



0.7599 



^^1=2.5125 ; ^2^2 = 1-9090 



0.6364 



1.8033 



1.6943 

9I2 = -1.6583 
g\i = 2.1994 



(69) 



The numerical values of parameters Qa can be obtained by means of Eqs (21), (22) and (24) 
where the distribution functions J-""*^ are constructed by using Eqs (61). From Eq. (21) we 



obtain an equation, depending on variable Bq,, for the particle density n^. 



na(x, 0c 



(ivJ-"''(e„,v,x) 



a = {e,i) 



(70) 



By injecting Eq. (70) into Eq. (22), we have 

/rfvvj^°0(e^,v,x) 



Ua(x, Oc 



/t/vJ^"O(0„,v,x) 



(71) 



This latter equation should then be combined with Eq. (24) and, for a given value of ASph 
we finally derive the equation for 6^ 

1 



-m„ / rfr I V - u„(x, e„) P j^°°(v, X, e„) = as, 



*Q:v ; ^-^aj 



, ^»., ^^aj 



Oh. 



a = e,« 



(72) 



with Ua is provided by Eq. (71 ). Eq. (72 ) can be solved (numerically) with respect to variable 
6q,; we find 

(Be, Qi) = (4.1760 X lO^^cm^sec-^ 2.2745 x lO^^cm^sec-^) (73) 



From Eq. (^70h, Eq. (71) and Eq. (23), we can now derive the profiles of the particle density. 



riaix), the average velocity, u^, and temperature, T^, respectively. The knowledge of the 
fields Ua, Ta etc. allows determining the profiles of the thermodynamic forces Z^ [j = (l,2,3) 
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for electrons, and j = (l,2) for ions]. Parameters P<^o and Aq are determined by Eqs (|9|. By 



substituting Eqs rtlGft into Eq. (19), we can easily check that distribution function, Eq. (19), 



tends to reduce to a pure gamma-process as P^f, — )■ ip{r), A — )■ (2 l-Bl)^"*^ and v^\ -^ { i.e., 
the pitch angle 6p is close to vr/2). This happens when the thermodynamic forces, X^ and 
X^, tend to vanish. Our aim is to derive the expressions of the steady-state electron and 
ion distribution functions for plasmas in the Onsager region, and fy, w, P^ and A close to 
0, EQ, ipir) and [2 |_B| (r,6') ] 7-*^ respectively (see also the footnote [36j). To this purpose, we 
shall work with the dimensionless variables P^ and A, which for the magnetic configuration 



given by Eq. (49), are defined as 



Boa? 
X = BoX = 

with p = 



Bn 



2 I B 



a^lloc \ B \ Jq q[p') 



w 
B 



(74) 



Bo{l- -—p cos 



We have 



00 



pa 



nii^O 



A" = A^ 



27r 



5n 



PaO 



P 



q{p' 
1 



hii^O 



- dp' = tpipao) 

-Paocosfc', 



aO 



a 



a 



[e,t} 
{e,i) 



(75) 



2\B\ 2 V Ro 
where (peo^^eo) and (piO)^io) are solutions of the equations (see the footnote [37]) 



Xlir, 



e(3) , -e(3) 



(76) 



0{e) \Xl{r,d)=g;i''>+gl'> = 0{e) 

X2(r,^) = gf^ = 0{e) ' | r,o G {r | x^l\r,9) = 0{e)} 

with the drift parameter e of the order e ~ 10^^ [12j. As an example of calculation, we choose 
the following values of the parameters: Bq = 3.45 Tesla, Rq = 2,96 m and a = 1, 25 m. In 



this case, the systems (76) admits the solutions 



(Peo,^eo) = (0.2376,1.0562) 

X] = X2 = 2.0367 X 10-^ ~ C(e) 



(Pio,M = (0.5079,2.8137) 
X} = 



(77) 



Fig. (12) illustrates the surfaces of the electron forces X^ and Xf. The intersection line 



corresponds to the values of {p,0) such that X^ = X^. Fig. (13) shows the curve p = p{6) 



where the ion thermodynamic force X} vanishes. The values of parameters P?g and Ag are 



easily estimated from Eqs ( 75 ) 
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X & X surface 





FIG. 12: The electron thermodynamic forces 
Xg and Xg. The intersection Hne corresponds to 
the values {p,0) for which X^ = X|. The first 



system of equations in Eqs ( 76 ) , is satisfied in the 
narrow region around to the values {peO,Oeo) = 
(0.2376,1.0562). 



FIG. 13: Locus of the points (r, 9) 
where Xf =0. As an example of calcu- 
lation, we choose the point {pio,9io) = 
(0.5079,2.8137) of this curve. 



(P|o, K) = (0-1651, 0.5246) ; (P^q, X^) = (0.5499, 0.3984) 
The parameters APfp, AAq, AAi can be obtained from the relations 



dxl 


dxl 


dp 


' de 


dX^ 


dxl 


dp 


' de 

PeO,9eO 


dX} 


dX} 


dp 


' de 

PiafiiO 



PeOfieO 



PeOfieO 



PiofiiO 



(0.0039, 1.5957 X 10" 



-0.0028,-7.3501 x 10 



(6.6596 X 10"^C(e^)) 



-5\ 



By the chain rule for derivatives, we find 



^1 ~ -3.9304 X 10"^AP| 
6 ^ 3.3571 X 10"^ AA' 



(78) 



(79) 



(80) 
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From which we obtain 



API = 182.278 



(API 



0{e 



AAg -> oo 



AAq -)■ oo 



AA^ = 581.268 
AAi = 286.236 



Hence, the electron and ion density distribution functions finally read 

'A-0.5246\2 



T-eO 

y- oc exp 
J"*" oc exp 



w-0.5105\2 

0.3096 ) 
10^ X ^-0.2780\2 



P-0.1651\2 
182.278 ) 



581.268 



\J\ 



A-0.3984\2 



\J\ 



(81) 



(82) 



0.1686 / V 286.236 

where the dimensionless variable w = w/v^f^^, with vthe computed at the center of the 



tokamak, has been introduced. Therefore, Eqs. (68) ensure that for collisional tokamak- 
plasmas in the Onsager region, the DDF J-'° coincides exactly with the reference DDF 
estimated by the neoclassical theory. In the coordinates w, P^ and A (and ?/'), the expression 



(1) 



of the DDF is given by Eq. (19). Notice that in this case Cg = (AAq -^ oo) and the 



.(0) 



presence of the parameter C2 (or of the parameter AAf ) is crucial. 



V. CONCLUSIONS 

Using statistical thermodynamics approach we have derived the general expression of the 
(density of) distribution function J-'*' for the case of a thermodynamic system out of equi- 
librium, subject to three thermodynamic forces. The local equilibrium is fixed by imposing 
the following conditions : 

i) The minimum entropy production condition on the two Prigogine's fiuctuations ai and 
"2; 

ii) The maximum entropy principle on the variable w, for ai = 0:2 = 0; 

ill) The scale invariance of the restrictions used in the maximization of the entropy; 

iv) A new mathematical ansatz, used in selecting a minimal number of restrictions and 
implicitly free parameters. 

From this ansatz results a singularity of the DDF that has immediate physical interpretation 
in terms of the intermittency in turbulent plasmas. 
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The derived DDF, J-'°, is more general than that currently used for fitting the numerical 
steady-state solution obtained by simulating ICRH plasmas and for describing various sce- 
narios of tokamak plasmas. The adopted procedure can be generahzed for systems subjected 
to an arbitrary number of thermodynamic forces. By kinetic theory, we have linked, and 
then fixed, the seven free parameters entering in J^ with the external energy sources and 
the (internal) entropy production source strength. To be more concrete, we have analyzed 
the case of, fully ionized, magnetically confined plasmas. 

A perspective of this work is the following. In 2009, a macroscopic theory for closure 
relations has been proposed for systems out of Onsager' s region [2S]- This theory is referred 
to as the thermodynamic field theory (TFT). The aim of the work is to determine the 
nonlinear flux-force relations that respect the thermodynamic theorems for systems far from 
equilibrium. Through the TFT it is possible to estimate the nonlinear DDFs for coUisional 
tokamak-plasmas when the nonlinear contributions cannot be neglected [29]. Also in this 
case, the crucial question is : When does our reference DDF identify with the, steady- 
state, nonlinear distribution function derived by the TFT ? The solution of this difficult 
problem will contribute to establish the link between the macroscopic approach (TFT) 
with a microscopic description based on statistical thermodynamics. Another problem to be 
solved is the possibility to improve the numerical fit by adding new free parameters according 
to the principles exposed in the Appendix. The work is ongoing. 
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Appendix A: The Physical Justification of the Gamma Distribution Function 

In order to justify the w dependence of V{w) we will use the principle of maximal entropy 
with suitable chosen restrictions obtained from scale invariance and a new mathematical 
ansatz. We denote by 5'[p(.)] the entropy functional of a given probability density function 
p{w) : 

S[p{.)] = - piw)\og{p{w))dw (Al) 

The maximum entropy principle will be applied by imposing a set of restrictions: 





oo 



p{w)dw = 1 (A2) 

wp{w)dw = E{w) = pi (A3) 

\og{w)p{w)dw = E{\og{w)) = V (A4) 





oo 



The naturalness of the restriction given by Eq. (A3) is clear: we observe that px = T/m 
where T is the temperature. Concerning the naturalness of the choice of restriction (A4), we 
remark that all of the previous restrictions belongs to a class that is invariant under scaling. 
More generally, for a restriction of the form E(w") = pa-, respectively for E(log(t(;)) = z/, the 
effect scaling w = kw' is E(u7'") = p'^ = pak~°', respectively for E(log(w')) = u' = u — log^k). 
The problem why the PDF obtained by the one of the possible simplest choice from the more 
general set of restrictions 

E («;"*=)= /ifc ; with /xq = 1; ao = 0; cti = 1; fc = 0, 1, 2, . . . , (A5) 

E(log(w)) = u (A6) 

deserves further study. A partial answer is given by simplicity and extremality reasons. 



Observe first that the restriction (A6) can be seen as a limiting case of (A5). Indeed, 



suppose that we have for some fixed k : a^ = e ^ 1 

E{w') = /ifc (A7) 



From (A2) and (A7) results 

E (^) = ^ (A8) 

When the support of the PDF p{w) is concentrated mainly on the domain when | log(w)| 
is not too large, then we can approximate: {w'^ — w^)/e = log{w) , so the Eq.(A8) is 
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reduced to Eq.(A4). So, the restriction of the type (A4) can be seen as the an extreme 



case. It is easily checked that applying the maximal entropy principle from X6 results 
p{w) X const w'^~^ for some 7 > 0. Notice that when 7 < 1, the singularity of our 

io->-0 

DDF is related to the intermittency shown by real physical DDFs [30], with time and 
ensemble average are provided by our DDF. On the other hand, it is easily checked that 



in the case of restrictions, A5 after applying the maximal entropy principle, we obtain 
p{w) X ci exp(— C2tf™) for some Ci,C2 > 0, where m = max{ai, ..., a„}. But m > 1 give 
a much faster decay compared to Maxwell distribution, so we obtain m < 1. Consequently, 



we use the minimal ansatz given by Eqs (A2)-(A4), for the selection of the restrictions. 



From this reason of extremality of the restriction Eq.(A4) and motivated by the minimality 



assumptions, we will explore the consequences of the minimal and extremal model given by 



Eqs (A2)-(A4). In order to find the distribution function p{w) that maximizes the entropy 



(Al), with restrictions (A2)-(A4) we use the Lagrange multiplier method. The result is the 



gamma distribution function, given by 



w/Qy exp (— w/0) 



(A9) 



where r(7) is the Euler gamma function. The relation between the parameters 7, G and 
Pi = T/m, V is given by 



E(w) = T jm = 7O 
E(log(ti;)) = z/ = ^(7) + log(e) 



(AlO) 
(All) 



where \&(7) is the digamma function. Due to the very special mathematical peculiarity in the 



ansatz Eqs (A2)-(A4), the resulting PDF is expected to have also some special properties. 



Indeed, the gamma distribution is infinitely divisible and stable : if the PDF of independent 
random variables Xi, . . . ,Xn is a gamma distribution with the same scale parameter and 
shape parameters 71, . . .7™, then PDF of the random variable X]fc=i^fc is again a gamma 
distribution with the same scale parameter and with shape parameter Ylt=i 7^- ^^ ^^^ 
convinced that the physical interpretation of the infinite divisibility property is a challenging 
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problem. 
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